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We revisit the scenario of the gravitational radiation recoil acquired by the final remnant of a 
black-hole-binary merger by studying a set of configurations that have components of the spin both 
aligned with the orbital angular momentum and in the orbital plane. We perform a series of 42 new 
full numerical simulations for equal-mass and equal-spin-magnitude binaries. We extend previous 
recoil fitting formulas to include nonlinear terms in the spins and successfully include both the new 
and known results. The new predicted maximum velocity approaches 5000km/s for spins partially 
aligned with the orbital angular momentum, which leads to an important increase of the probabilities 
of large recoils in generic astrophysical mergers. We find non-negligible probabilities for recoils of 
several thousand km/s from accretion-aligned binaries. 
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Introduction: With the breakthroughs of 2005 in the nu- 
merical techniques to evolve black-hole binaries (BHBs) 
[TH3] , Numerical Relativity (NR) became a very impor- 
tant tool to explore highly-dynamical and nonlinear pre- 
dictions of General Relativity. In the last few years we 
have gained notable insight into the modeling of gravi- 
tational radiation to assist laser interferometer detectors 
0]. There are also numerous examples of explorations 
in the realm of Mathematical Relativity and in the as- 
trophysical scenarios for supermassive black-hole mergers 
and retention of black holes in galaxies and globular clus- 
ters 01]. 

Some of the most striking recent discoveries are related 
to effects due to the intrinsic spin of the individual BHs 
during the final merger stage of BHBs. In general, the 
spins and orbital plane precess during the inspiral and 
the spin direction of the remnant is misaligned with the 
individual BH spins [7 prior to merger (spin flips) . Spins 
can also have a dramatic effect on the inspiral rate. When 
the BH spins are partially aligned with the orbital angu- 
lar momentum, the merger is delayed, while when they 
are antialigned, the merger happens much more quickly. 
This "hangup" effect [8] is due to an unexpectedly strong 
spin-orbit coupling. Perhaps even more surprisingly are 
the very large recoils [5] acquired by the remnant of the 
merger of comparable-masses, highly-spinning, BH in the 
"superkick" configuration, where the spins lie along the 
orbital plane (equal magnitude, but opposite in direc- 
tion) [51 |H1 HO] ■ Initial studies, which indicated that these 
BHBs recoil at up to ~ 4000 kms" 1 IP, prompted 
astronomers to search for possible recoil candidates. To 
date, a few interesting cases of galaxies with cores dis- 
playing differential radial velocities of several thousands 
km/s [T2TTT4] have been found. More systematic recent 
studies produced tens of potential candidates [UJ [TB] . 



The discovery of large recoils also triggered theoretical 
statistical studies of BHB dry mergers [T7] and Monte 
Carlo simulations of wet premergers to study the effect 
of accretion and resonances on spin distributions |18l - 
[20] . Accretion tends to align spins with the orbital angu- 
lar momentum [211 122] ; resulting in a notable reduction 
in the probabilities of observing large recoils, either di- 
rectly, or through their influence on the galactic cores, 
for small and medium sized galaxies (SJ. For BHs with 
masses larger than 10 s M©, alignment by accretion is less 
effective. In this Letter we revisit the scenarios for the 
generation of recoils by studying a set of configurations 
that combine two of the largest spin effects observed in 
BHBs, the hangup effect and superkicks. The combined 
effect appears to be a dramatic increase in the probability 
distribution for large recoils, as we will show below. 

Full Numerical Simulations: We evolved a set of 42 
equal-mass, spinning, quasicircular configurations using 
the LazEv |23] implementation of the moving puncture 
formalism [3] , with the conformal factor W — ^fx = 
cxp(— 2(f)) suggested by [24 as a dynamical variable. For 
the runs presented here we use centered, eighth-order fi- 
nite differencing in space |25] and a fourth-order Runge- 
Kutta time integrator. The LazEv code used the Cac- 
tus/EinsteinToolkit [26J [27] numerical infrastructure 
along with the Carpet [55] mesh refinement driver. We 
use the TwoPunctur.es [55] thorn to calculate the ini- 
tial data. We use AHFinderDirect [SU] to locate ap- 
parent horizons. We measure the magnitude of the hori- 
zon spin using the Isolated Horizon algorithm detailed 
in [31]. 

Our configurations have the property that the in- 
plane components of the spins of the two BHs have the 
same magnitude, but opposite signs, while the out-of- 
plane components have the same magnitude and sign (see 
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FIG. 1: Black-hole-binary configuration for hangup recoils. 

Fig. fTl). They thus combine the hang-up [8] and super- 
kick [51 [3] effects. Additionally, the orbital plane does not 
precess, but rather moves up and down along the direc- 
tion of orbital angular momentum as the binary evolves. 

We performed a set of 30 simulations with individual 
BH spins of magnitude a = l/v2 and 12 simulations 
with BH spin magnitudes of a = 0.91, where a is the 
normalized spin of the BH. The a = l/y/2 conn gurations 
were split into five sets of 6, where the runs in each indi- 
vidual set had the same initial angle 9 between the spin 
direction and orbital angular momentum direction (here 
we chose 8 = 22.5°, 45°, 60°, 120°, 135°). In each set 
with a given 8, we chose the initial orientation </>j be- 
tween the in-plane spin and linear momentum to be 0°, 
30°, 90°, 130°, 210°, and 315°. For the a = 0.91 runs, 
we used the same initial 6 <pi configurations for 8 = 60° 
and 8 = 15°. We combine these results with the simu- 
lations of [IT] (which have 8 = 90°) in order to perform 
our analysis below. 

We set up the initial separations, such that each binary 
completed 5-6 orbits, prior to merger (to reduce eccen- 
tricity). The initial separations varied between 10.16M 
and 8.2M (depending on the magnitude of the hangup 
effect). 

Results and Analysis: In a previous study we found 
that the superkick recoil (where the two BHs have equal 
mass, equal intrinsic spin magnitudes a, and spins lying 
in the orbital plane in opposite directions) has the follow- 
ing dependence on spin a and orientation cf> (the angle 
between the in-plane spin vector and the infall direction 
near merger), 

V = Vicos(<t>-<l)i) + V 3 co8(3<t>-34> 3 ), 
Vi = Vi S a + V h3 a 3 , 

V 3 = V 3A a + V 3 , 3 a 3 , (1) 

where V h3 = (-15.46 ± 2.66) kms" 1 , V 3 .i = (15.65 ± 
3.01) kin's -1 , and V 3 , 3 = (105.90 ± 4.50) kms" 1 , while 
Vi.i = (3681.77 ± 2.66) kms -1 . From that study, it was 
clear that in the superkick configuration, the dominant 
contribution, even at large a, is linear in a and propor- 
tional to cos(c/)). Note that because of the small con- 
tributions of V 3 and Vi.,3, we neglect these terms in the 
statistical studies below (where we take a uniform distri- 
bution in <f> — <pi). 

Our initial motivation for the current study was to 
determine if the hangup effect [5], which amplifies the 



amount of radiation emitted by the BHB, also affects the 
maximum recoil. To this end, we looked at configurations 
that combined both the hangup and superkick effects. 
Based on the superkick formula ([lj, we expected that 
the recoil would have the form 

Vi = Vi^a sin 6 + Aa 2 sin 8 cos 8 + 
Ba 3 sin 8 cos 2 8 + Ca 4 sin 9 cos 3 9, (2) 

where V\ is the component of the recoil proportional to 
cos0, Vj. i arises from the superkick formula, and the 
remaining terms are proportional to linear, quadratic, 
and higher orders in S z /m 2 = a cos 8 (the spin com- 
ponent in the direction of the orbital angular momen- 
tum). Here, we do not consider terms higher-order in 
the in-plane component of A oc d?2 — qa\ denoted by A 1 - 
(A 1 - oc a sin 9 here), where q = m\/rri2 is the mass ra- 
tio, because our previous studies showed that these terms 
were small at 9 = 90°. A fit to this ansatz ^ showed 
that the coefficients V ls = (3677.76± 15.17) kms" 1 , A = 
(2481.21 ± 67.09) kms" 1 , B = (1792.45 ± 92.98) kms" 1 , 
C = (1506.52 ± 286.61) kms -1 converge very slowly and 
have relatively large uncertainties. In addition, we pro- 
pose the modification 

„ . / 1 + Ea cos 9 \ 

V x = Dasind , 3 

V 1 + Fa cos 8 J w 

which can be thought of as a resummation of Eq. |2]) 
with an additional term EacosS, and fit to D, E, F 
(where we used the prediction of [TT] to model the V\ for 
9 = 90°) and find D = (3684.73 ± 5.67) kms" 1 , E = 
0.0705 ±0.0127, and F = -0.6238 ±0.0098. Note that E 
is approximately 1/10 of F, indicating that corrections 
to this formula converge quickly. The two formulas ^ 
and ([3]) give very similar results for a broad range of a. 
We then use Eq. ^ to predict the recoil for higher spin 
a = 0.91 and test this formula for three angles 8 = 90°, 
8 = 60°, and 8 = 15°, with very good agreement (see 
Fig. [2]). In actuality, both Eq. Q and Eq. Q provide 
accurate predictions for our measured recoils at a = 0.91. 
The results are startling. The recoil is not maximized at 
8 = 90° , as was previously assumed based on linear spin- 
orbit PN expressions, but rather at smaller angles (see 
Table [T]) . Additionally, the maximum recoil is closer to 
5000km s _1 . Note that we measured recoils from actual 
simulations as large as 4171km s" 1 for a = 0.91, which is 
larger than the previously predicted maximum possible 
recoil of 3681km s _1 for quasicircular binaries. This can 
have profound astrophysical implications because partial 
alignment of the binary (e.g. due to accretion), rather 
than inhibiting large recoils, can actually amplify them, 
leading to much larger probabilities for observing high 
recoils. 

Using the same post-Newtonian analysis [52] as in 33! , 
we can extend formulas ([2]) and (J3J) to less symmetric con- 
figurations by replacing a sin 8 by [ajr — qa~^]/(l + q) and 
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TABLE I: The angle 8 that gives the largest recoil along with 
the magnitude of this recoil for different values of the indi- 
vidual BH spin a. The columns to the left are from Eq. (3), 
while the columns to the right are for Eq. 



a 


Pmax 


Knax 


^raax 


Knax 


0.1 


86.02° 


369.36 


86.13° 


368.62 


0.5 


70.16° 


1961.38 


69.99° 


1955.51 


1/V2 


61.90° 


2968.71 


61.33° 


2967.85 


0.91 


53.55° 


4224.93 


53.92° 


4231.93 


1 


49.67° 


4925.94 


51.22° 


4915.22 



TABLE II: Probability that the recoil velocity will be in a 
given range P, and the probability that the recoil will be a in 
a given range along the line of sight P ODS (to relate to possi- 
ble redshift measurements in galaxies) for the new prediction 
(LEFT) and the old prediction (RIGHT) 



range 


P 


-fobs 


P old 


P obs old 


0-500 


79.027% 


92.641% 


94.888% 


98.914% 


500-1000 


15.399% 


6.177% 


4.921% 


1.067% 


1000-2000 


5.384% 


1.164% 


0.191% 


0.019% 


2000-3000 


0.189% 


0.018% 








3000-4000 


0.001% 


0.0001% 









FIG. 2: A fit of the recoil (Vi) to the form Eq. g for the 
a — l/v2 configurations, and predictions (based on this fit- 
ting) for the a — 0.91 recoils. Note how well the a = 0.91 
curve matches the three measured values. For reference, 
curves corresponding to the original empirical formula pre- 
diction (which only had terms linear in A) for a = 1/ v2 and 
the new formula for a — 1 are also included. Note the skew 
in the velocity profile compared to the linear predictions. 

5000 





a=0.91 (Nonlinear) 




N. ct=0.707 (Nonlinear) 




' >. ct=1 (Nonlinear) 




\. \^ a=0.707 (Linear) 







a cos 9 by 2[a| 



q 2 a\ 



]/(l + 9) 2 



Importantly. 

„_L 



we are 

assuming that terms proportional to [a^ — qct\] n (for 
n > 1) are negligible. This can be verified by confirming 
that formulas (|2| and ^ are accurate for all 9 and a (a 
subject of our ongoing analysis that will be reported in 
a forthcoming paper). We emphasize that the proposed 
extension is an ansatz, that while reasonable as a start- 
ing point for the modeling, needs to be thoroughly tested 
and refined. In this letter we will use Eq. ^ when gener- 
alizing to unequal masses and arbitrary spin orientations. 
A generalization of the resummation ([3| to the generic 
mass ratio q case is nontrivial, and will be discussed in an 
upcoming paper. Our ansatz for the generic recoil (for 
brevity, we only display the term in the direction of the 
orbital angular momentum proportional to Aj_) then is 



16- 



(1 + 9) 1 ' 

C0S(<^>A - <Al), 



qai\ 



V 1 . 1 +AS + BS 2 + CS- : 



(4) 



where S — 2(a| 



h q 2 af)/(l + q) 2 , <j>A is angle between 
and the infall direction (evaluated at 
a fiducial point around merger), <j>\ is a constant, and 



1_L oc oti 



qa^ 



FIG. 3: The recoil probability distribution using the new and 
old empirical formulas for the recoil, starting from a distri- 
bution of BHB configurations consistent with recent models 
for accreting binaries [37] . The new formula predicts a signif- 
icantly larger probability for high recoils. 

P(v) 




r) = q/(l + q) 2 is the symmetric mass ratio. See Ref. [T7] 
for the remaining components of the recoil velocity. 
Astrophysical Implications: To test the effect of our 
new empirical formula on the predicted recoil rates, we 
consider a model distribution of BHBs with mass ra- 
tio distribution P(q) cx q~ a 3 (l — q) [341436] . spin dis- 
tribution P(a) oc (1 — a)^ 1 ' cr 0-1 ' , with parameters 
a = 4.8808 and b = 1.72879 (which models the distribu- 
tion kindly provided by M. Volonteri [3 7J resulting from 
spin- up effects due to accretion), spin-direction distribu- 
tion P(9) oc (1 - 9) {b - l) 9 {a - 1 \ with parameters a = 2.5 
and 6 = 7 (which approximates the spin-direction distri- 
bution provided by M. Volonteri [37j . which is a simpli- 
fied model for the spin distribution for hot accretion 22 ) , 
and P{4>) uniform in < < 2tt. We model 10 7 BHBs 
consistent with these distributions and measure the pre- 
dicted recoil using Eq. Q, as well as our original em- 
pirical formula, and compare the probabilities for ob- 
serving a recoil in a given velocity range. Our results 
are summarized in Fig. [3] and Table [TTJ The dramatic ef- 
fect of the new recoil prediction is apparent. Recoils in 
the range 1000 kms -1 — 3000km s _1 are significantly en- 
hanced, leading to a realistic chance for observing large 
recoils. Details on how the statistical studies were per- 
formed can be found in [17j . 

Conclusions and Discussion: We revisited the scenario 
for the generation of large gravitational radiation recoils 
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acquired by the remnant BH after the merger of BHBs 
and found that configurations with spins partially aligned 
with the orbital angular momentum produce larger re- 
coils (up to 1200 km/s more) than those with spins lying 
in the orbital plane (aka superkicks) . The new configura- 
tion maximizes the total momentum radiated by optimiz- 
ing over the competing requirements of maximizing the 
total power radiated (which occurs for the hangup config- 
uration) and the skew in the power distribution (which 
is maximized in the superkick configuration). Our re- 
sults imply a nonlinear coupling among components of 
the spins that can be expressed in the simple form given 
by Eq. Q. Based on these new terms in the empirical 
formula for recoils, we recalculate the probabilities for 
large recoils to occur in astrophysical scenarios of BHB 
encounters when accretion effects are included. At small 
angles of the spins with the orbital angular momentum, 
the new term magnifies recoil velocities by up to a fac- 
tor 2.8 with respect to the previous formula, and we find 
non-negligible probabilities of observing black holes re- 
coiling at several thousand km/s, as reported in Table [TT] 
Our results indicate that there is a need for additional 
theoretical searches for large recoils in other regions of 
the parameter space, and lend support for additional ob- 
servational searches for high-velocity black holes. 
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